# TODO: Add comment
# 
# Author: yaping
###############################################################################
##parameters: working_dir, file_one, file_two, NDR_p_value_threshold, annotation_file, overlap_or_not_with_annotation_location, output_bed_file_in_different_combination_with_gene_names
##each column's content:  chr	start	end	aveMethyWind	gchNumWind	gchDepthWind	gchCTdepthWind	gchDotWind	sigValue	aveMethyWindLinker	wcgMethyWind	wcgNumWind


for (e in commandArgs(TRUE)) {
	ta = strsplit(e,"=",fixed=TRUE)
	if(! is.na(ta[[1]][2])) {
		if(ta[[1]][1] == "wd"){
			wd<-ta[[1]][2]  ## directory to make plot file
			
		}
		if(ta[[1]][1] == "file1"){
			file1<-ta[[1]][2]
		}
		if(ta[[1]][1] == "sample1"){
			sample1<-ta[[1]][2]
		}
		
		if(ta[[1]][1] == "pValue"){
			pValue<-ta[[1]][2]
		}
		if(ta[[1]][1] == "anno"){
			anno<-ta[[1]][2]
		}
		if(ta[[1]][1] == "annoName"){
			annoName<-ta[[1]][2]
		}
		if(ta[[1]][1] == "annoFormat"){
			annoFormat<-ta[[1]][2]
		}
		if(ta[[1]][1] == "annoRegionExtend"){
			annoRegionExtend<-ta[[1]][2]
		}
		if(ta[[1]][1] == "overlap"){
			overlap<-as.numeric(ta[[1]][2])  ##0: no overlap; 1: overlap
		}
		if(ta[[1]][1] == "outputBed"){
			outputBed<-as.numeric(ta[[1]][2])  ## 1: output; others: no
		}
		if(ta[[1]][1] == "useWcg"){
			useWcg<-as.numeric(ta[[1]][2])  ## 1: use cg methylation value; others: no
		}
		if(ta[[1]][1] == "dataPoint"){
			dataPoint<-as.numeric(ta[[1]][2])  ## 1: use cg methylation value; others: no
		}
		if(ta[[1]][1] == "gchNum"){
			gchNum<-as.numeric(ta[[1]][2])  ## 1: use cg methylation value; others: no
		}
	}
}
setwd(wd)
library(Ringo)

Obj1<-read.table(file1,sep="\t",header=F)


Obj1<-Obj1[!is.na(Obj1[,9]) & !is.na(Obj1[,4]) & Obj1[,8] > dataPoint & Obj1[,5] > gchNum,]


Obj1 <-cbind(Obj1, p.adjust(Obj1[,9],method="BH"))




Obj1Data<-Obj1



Obj1NDRInObj1<-Obj1Data[Obj1Data[,13]< pValue,]






if(!is.null(anno)){
	
	annoFile<-read.table(anno,sep="\t",header=F)
	#annoFile<-annoFile[annoFile[,6]>=50,]
	#annoFile<-annoFile[annoFile[,8]>5,]
	annoLoc<-data.frame(chr= annoFile[,1], start= annoFile[,2]-as.numeric(annoRegionExtend), end=annoFile[,3]+as.numeric(annoRegionExtend))
	if(annoFormat == "gtf"){
		annoLoc<-data.frame(chr= annoFile[,1], start= annoFile[,4], end=annoFile[,5])
	}
	callableWindowObj1Data=Obj1Data
	
	callableWindowLocObj1<-data.frame(chr= callableWindowObj1Data[,1], start= callableWindowObj1Data[,2], end=callableWindowObj1Data[,3])
	
	
	
	Obj1DataAnnoInfo<-array()
	AnnoDataAnnoInfo<-array()
	
	for(i in c(1:22)){
		chrName=paste("chr",i,sep="")
		tmpObj1<-callableWindowObj1Data[callableWindowObj1Data[,1] %in% chrName,]
		
		
		tmpAnno<-annoFile[annoFile[,1] %in% chrName,]
		
		tmpObj1NDR<-tmpObj1[tmpObj1[,13]<pValue,]
		
		Obj1LocData<-data.frame(chr= tmpObj1[,1], start= tmpObj1[,2], end= tmpObj1[,3])
		
		
		annoLocTmp<-data.frame(chr= tmpAnno[,1], start= tmpAnno[,2], end= tmpAnno[,3])
		
		tmpObj1NDRLocData<-data.frame(chr= tmpObj1NDR[,1], start= tmpObj1NDR[,2], end= tmpObj1NDR[,3])
		
		if(annoFormat == "gtf"){
			annoLocTmp<-data.frame(chr= tmpAnno[,1], start= tmpAnno[,4], end=tmpAnno[,5])
		}
		
		Obj1AnnoOverlapped<-as.matrix(regionOverlap(Obj1LocData, annoLocTmp))
		Obj1NDRAnnoOverlapped<-as.matrix(regionOverlap(annoLocTmp, tmpObj1NDRLocData))
		
		Obj1DataAnnoInfo<-rbind(Obj1DataAnnoInfo,cbind(tmpObj1, rowSums(Obj1AnnoOverlapped)))
		AnnoDataAnnoInfo<-rbind(AnnoDataAnnoInfo,cbind(tmpAnno, rowSums(t(Obj1AnnoOverlapped)), rowSums(Obj1NDRAnnoOverlapped)))
		
	}
	Obj1DataAnnoInfo<-Obj1DataAnnoInfo[2:length(Obj1DataAnnoInfo[,1]),]
	AnnoDataAnnoInfo<-AnnoDataAnnoInfo[2:length(AnnoDataAnnoInfo[,1]),]
	
	Obj1Anno<-Obj1DataAnnoInfo[Obj1DataAnnoInfo[,14]!=0,]
	
	
	Obj1NDRannoInObj1<-Obj1DataAnnoInfo[Obj1DataAnnoInfo[,14]!=0 & Obj1DataAnnoInfo[,13]<pValue,]
	
	
	
	if(overlap==0){
		
		
		Obj1Anno<-Obj1DataAnnoInfo[Obj1DataAnnoInfo[,14]==0,]
		
		Obj1NDRannoInObj1<-Obj1DataAnnoInfo[Obj1DataAnnoInfo[,14]==0 & Obj1DataAnnoInfo[,13]<pValue,]
		
		
		
	}
	
	
	
}


cat("Total callable window of Obj1:")
cat(dim(Obj1)[1])
cat("\n")

cat("Total NDR of Obj1:")
cat(dim(Obj1NDRInObj1)[1])
cat("\n")

cat("NDR percentage in callable window of Obj1:")
cat(dim(Obj1NDRInObj1)[1]/dim(Obj1)[1])
cat("\n")

if(!is.null(anno)){
	cat("Total callable window of Obj1 in annotation file:")
	cat(dim(Obj1Anno)[1])
	cat("\n")
	
	cat("Percentage of callable window of Obj1 in annotation file:")
	cat(dim(Obj1Anno)[1]/dim(Obj1)[1])
	cat("\n")
	
	cat("Total NDR of Obj1 in annotation file:")
	cat(dim(Obj1NDRannoInObj1)[1])
	cat("\n")
	
	cat("NDR percentage in callable window of Obj1 in annotation file:")
	cat(dim(Obj1NDRannoInObj1)[1]/dim(Obj1Anno)[1])
	cat("\n")
	
	cat("NDR percentage in NDR of Obj1:")
	cat(dim(Obj1NDRannoInObj1)[1]/dim(Obj1NDRInObj1)[1])
	cat("\n")
	
	cat("(anno+NDR)/anno:")
	cat((dim(Obj1Anno)[1] + dim(Obj1NDRInObj1)[1])/dim(Obj1Anno)[1])
	cat("\n")
	
	cat("(anno+NDR)/NDR:")
	cat((dim(Obj1Anno)[1] + dim(Obj1NDRInObj1)[1])/dim(Obj1NDRInObj1)[1])
	cat("\n")
	
}

if(outputBed == 1){
	filename=paste(sample1,annoName,"pValue",pValue,"DP",dataPoint,"GCH",gchNum, "summary.txt",sep=".")
	summary<-cbind(dim(Obj1)[1],dim(Obj1NDRInObj1)[1],dim(Obj1Anno)[1],dim(Obj1Anno)[1]/dim(Obj1)[1],dim(Obj1NDRannoInObj1)[1],dim(Obj1NDRannoInObj1)[1]/dim(Obj1Anno)[1],(dim(Obj1Anno)[1] + dim(Obj1NDRInObj1)[1])/dim(Obj1Anno)[1],(dim(Obj1Anno)[1] + dim(Obj1NDRInObj1)[1])/dim(Obj1NDRInObj1)[1])
	colnames(summary)<-c("Obj1_callableWindow","Obj1_NDRs","Anno_in_Obj1_callableWindow","Perc_Anno_in_Obj1_callableWindow","Obj1_NDR_overlap_anno","Perc_Obj1_NDR_overlap_anno","(anno+NDR)/anno","(anno+NDR)/NDR")
	write.table(summary,filename,quote=F,sep="\t",row.names=F,col.names=T)
}


